---
title: "Uso do R nas Pesquisas"
author: "Patrícia Pain"
date: "25/04"
output:
  html_document: default
  word_document: default
  pdf_document: default
---

# Por que o R?

## Uma ex-"StataUser"

2018 (TCC) - 2022 (Dissertação)

[![Fonte: https://www.stata.com/stata12/interface/](images/clipboard-2789924045.png)](https://www.stata.com/stata12/interface/)

## Por que o R?

![](images/clipboard-3229384149.png){width="100"}

![](images/clipboard-2686762934.png){width="250"}

![](images/clipboard-2585036236.png){width="200"}

![](images/clipboard-2152667611.png){width="100"}

[![\<https://danmrc.github.io/R-para-Economistas/index.html\>](images/clipboard-2833022016.png)](https://danmrc.github.io/R-para-Economistas/index.html)

-   Migrar de uma ferramenta para outra custa tempo e esforço!

## Quanto tempo é muito tempo?

-   Disciplinas na UFES

-   Cursos

-   **PRÁTICA**

![](images/clipboard-3318222554.png)

![](images/clipboard-774168444.png){width="315"}

![](images/clipboard-2760514626.png){width="330"}

![](images/clipboard-1653830451.png)

![](images/clipboard-944534974.png){width="255"}

# Começando...

## Instalar e Carregar Pacotes

### Forma mais rápida

1º: criar um objeto com todos os pacotes que irá utilizar

```{r pacotes, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
pacotes <- c("AER","base", "basictabler", "bibtex","BiocManager", "bookdown","caper", "car","caret", "conflicted","correlation","corrplot","cowplot","DescTools","distill","dplyr","factoextra","FactoMineR", "faraway","fastDummies","flextable","foreign","gdata", "ggrepel","ggplot2","ggpubr","graphics", "grid","gridExtra","gtsummary","Hmisc","httr2","jsmodule", "jtools","knitr","kableExtra","knitLatex","lmtest", "lubridate","magick","margins", "marginaleffects", "MASS", "MatchIt","mfx","mgcv","minqa","modelr" ,"mgcv","nnet","nortest","OddsPlotty","papaja", "pandoc","palmerpenguins","performance", "pglm","plm","plotly","plotly", "pROC","pscl","psych","rddtools","readr", "regclass", "readxl","RefManageR", "remotes", "reshape2",  "report","reshape2","rgl","rlang","rmarkdown","Rmisc","ROCR","RSelenium", "scales","sjlabelled", "stargazer","stats","stringr", "stringi","texreg","tidyr","tidyverse", "tinytex","tseries","truncnorm", "visreg","viridis", "xtable","wesanderson", "xfun", "writexl")
```

2º: a função abaixo irá instalar (caso já não esteja instalado) e carregar os pacotes do objeto "pacotes" criado no código anterior

```{r instalando, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
options(repos = "https://cran.rstudio.com/")
if(sum(as.numeric(!pacotes %in% installed.packages())) != 0){
  instalador <- pacotes[!pacotes %in% installed.packages()]
  for(i in 1:length(instalador)) {
    install.packages(instalador, dependencies = T)
    break()}
  sapply(pacotes, require, character = T)
} else {
  sapply(pacotes, require, character = T)
}
```

## Definir diretório

Session \> Set Working Directory \> Choose directory...

```{r definindo diretorio, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
setwd("C:/Users/Usuário/OneDrive/Oficina UFSM 2024")
```

## Carregar a base

```{r}
load("C:/Users/Usuário/OneDrive/Oficina UFSM 2024/envir.RData")
```

# Tratamento de dados

Recomendo o vídeo em: <https://www.youtube.com/watch?v=AoSKH0x8U3Q&t=2553s>

*O script está em um link do Harvard Dataverse na descrição do vídeo.*

# Estatística Descritiva

```{r descritiva, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
summary(Base_GECAT)
```

```{r descritiva única variavel, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
summary(Base_GECAT$SIZE_w1)
```

### Por grupos

Comparando as variáveis transformadas de tamanho da empresa entre países desenvolvidos e emergentes (usando across)

```{r comparando, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
EstatisticaDescritivaSIZE <- Base_GECAT %>%
  group_by(Emerg) %>%
  dplyr::summarise(across(starts_with("SIZE_w1"), list(mean = mean, sd = sd), na.rm = TRUE)) %>%
  ungroup()

EstatisticaDescritivaSIZE
```

## tbl_summary

<https://www.danieldsjoberg.com/gtsummary/reference/tbl_summary.html>

-   For **categorical variables** the following statistics are available to display.

    -   {n} frequency

    -   {N} denominator, or cohort size

    -   {p} formatted percentage

-   For **continuous variables** the following statistics are available to display.

    -   {median} median

    -   {mean} mean

    -   {sd} standard deviation

    -   {var} variance

    -   {min} minimum

    -   {max} maximum

    -   {sum} sum

    -   {p##} any integer percentile, where \## is an integer from 0 to 100

    -   {foo} any function of the form foo(x) is accepted where x is a numeric vector

```{r sumarizando2, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
#pacotes necessários para esse comando: tidyverse, gtsummary, flextable

Base_GECAT %>%
  dplyr::select(Emerg, SIZE_w1, BTM_w1, COV, Country, LCS, AUD) %>%
  gtsummary::tbl_summary(statistic = list(all_continuous() ~ "{median}, {mean} ({sd}) {p25}, {p75}",
                                          all_categorical() ~ "{n} ({p}%)"
                                          ),
                         digits = all_continuous() ~ 2,
                         missing = "no", #"ifany" = mostra missing se tiver missing, "always" = quant. de missing na variável
                         missing_text = "Missing",
                         ) %>%
  add_n() %>% # Adicionando coluna de não missing
  modify_header(label ~ "**Variáveis**") %>%
  modify_caption("Tabela X. Estatística Descritiva") %>%
  italicize_levels() %>%
  as_flex_table() 
```

### Em uma janela temporal

```{r sumarizando3, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
Base_GECAT %>%
  dplyr::select(Year, Emerg, SIZE_w1, BTM_w1, COV, Country, LCS, AUD) %>%
  dplyr::filter(Year > 2018 & Year< 2022) %>%
  gtsummary::tbl_summary(statistic = list(all_continuous() ~ "{median}, {mean} ({sd}) {p25}, {p75}",
                                          all_categorical() ~ "{n} ({p}%)"
                                          ),
                         digits = all_continuous() ~ 2,
                         missing = "no", #ifany = mostra missing se tiver missing, always = quant. de missing na variável
                         missing_text = "Missing",
                         ) %>%
  add_n() %>% # Adicionando coluna de não missing
  modify_header(label ~ "**Variáveis**") %>%
  modify_caption("Tabela X. Estatística Descritiva") %>%
  italicize_levels() %>%
  as_flex_table() 
```

Renomeando as variáveis no output

```{r renomeando_output, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
Base_GECAT %>%
  dplyr::select(Emerg, SIZE_w1, BTM_w1, COV, Country, LCS, AUD) %>%
  gtsummary::tbl_summary(statistic = list(all_continuous() ~ "{median}, {mean} ({sd}) {p25}, {p75}",
                                          all_categorical() ~ "{n} ({p}%)"
                                          ),
                         digits = all_continuous() ~ 2,
                         missing = "no", #ifany = mostra missing se tiver missing, always = quant. de missing na variável
                         missing_text = "Missing",
                         label = list(Emerg ~ "Emergentes", 
                                      SIZE_w1 ~ "Tamanho",
                                      BTM_w1 ~ "BTM",
                                      COV ~ "Covid",
                                      Country ~"País",
                                      LCS ~ "Ciclo de Vida", 
                                      AUD ~ "Auditor")
                         ) %>%
  add_n() %>% # Adicionando coluna de não missing
  modify_header(label ~ "**Variáveis**") %>%
  modify_caption("Tabela X. Estatística Descritiva") %>%
  italicize_levels() %>%
  as_flex_table() 
```

### Por grupo de tratamento(Emergentes) e controle(Desenvolvidos) - 2 grupos

Gerando tabela

```{r tabela, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
Base_GECAT %>%
  dplyr::select(Emerg, SIZE_w1, BTM_w1, COV, Country, LCS, AUD) %>%
  dplyr::mutate(Emerg = factor(Emerg, labels = c(0, 1))) %>%
  gtsummary::tbl_summary(by = Emerg,
                         statistic = list(all_continuous() ~ "{median}, {mean} ({sd}) {p25}, {p75}",
                                          all_categorical() ~ "{n} ({p}%)"),
                         digits = all_continuous() ~ 2,
                         missing = "no", #ifany = mostra missing se tiver missing, always = quant. de missing na variável
                         missing_text = "Missing",
                         label = list(Emerg ~ "Emergentes", 
                                      SIZE_w1 ~ "Tamanho",
                                      BTM_w1 ~ "BTM",
                                      COV ~ "Covid",
                                      Country ~"País",
                                      LCS ~ "Ciclo de Vida", 
                                      AUD ~ "Auditor")
                                      ) %>%
  add_n() %>% # Adicionando coluna de não missing
  add_p(list(all_continuous() ~ "t.test",
             all_categorical() ~ "kruskal.test")) %>%
  modify_header(label ~ "**Variáveis**") %>%
  modify_caption("Tabela X. Estatística Descritiva") %>%
  italicize_levels() %>%
  bold_p() %>%
  as_flex_table() 
```

Exportando em word

```{r exportando, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
Base_GECAT %>%
  dplyr::select(Emerg, SIZE_w1, BTM_w1, COV, Country, LCS, AUD) %>%
  dplyr::mutate(Emerg = factor(Emerg, labels = c(0, 1))) %>%
  gtsummary::tbl_summary(by = Emerg,
                         statistic = list(all_continuous() ~ "{median}, {mean} ({sd}) {p25}, {p75}",
                                          all_categorical() ~ "{n} ({p}%)"),
                         digits = all_continuous() ~ 2,
                         missing = "no", #ifany = mostra missing se tiver missing, always = quant. de missing na variável
                         missing_text = "Missing",
                         label = list(Emerg ~ "Emergentes", 
                                      SIZE_w1 ~ "Tamanho",
                                      BTM_w1 ~ "BTM",
                                      COV ~ "Covid",
                                      Country ~"País",
                                      LCS ~ "Ciclo de Vida", 
                                      AUD ~ "Auditor")
                                      ) %>%
  add_n() %>% # Adicionando coluna de não missing
  add_p(list(all_continuous() ~ "t.test",
             all_categorical() ~ "kruskal.test")) %>%
  modify_header(label ~ "**Variáveis**") %>%
  modify_caption("Tabela X. Estatística Descritiva") %>%
  italicize_levels() %>%
  bold_p() %>%
  as_flex_table() %>%
  save_as_docx("Tab2", path = "descritiva2grupos.docx")
```

### Por País - 3 grupos

```{r exportando2, echo=TRUE, message=FALSE, warning=FALSE, paged.print=FALSE}
Base_GECAT %>%
  dplyr::select(Emerg, SIZE_w1, BTM_w1, COV, Country, LCS, AUD) %>%
  dplyr::mutate(Country = factor(Country, labels = c("Brazil", "Argentina", "Canada"))) %>%
  gtsummary::tbl_summary(by = Country,
                         statistic = list(all_continuous() ~ "{median}, {mean} ({sd}) {p25}, {p75}",
                                          all_categorical() ~ "{n} ({p}%)"),
                         digits = all_continuous() ~ 2,
                         missing = "no",
                         missing_text = "Missing",
                         label = list(Emerg ~ "Emergentes", 
                                      SIZE_w1 ~ "Tamanho",
                                      BTM_w1 ~ "BTM",
                                      COV ~ "Covid",
                                      Country ~"País",
                                      LCS ~ "Ciclo de Vida", 
                                      AUD ~ "Auditor")
                                      ) %>%
  add_n() %>% # Adicionando coluna de não missing
  add_p(list(all_continuous() ~ "t.test",
             all_categorical() ~ "kruskal.test")) %>%
  modify_header(label ~ "**Variáveis**") %>%
  modify_caption("Tabela X. Estatística Descritiva") %>%
  italicize_levels() %>%
  bold_p() %>%
  as_flex_table() 
#%>%
 # save_as_docx("Tab2", path = "descritiva3grupos.docx")
```

# Correlação

## Montar base para correlação só com as variáveis que quero

Todas as variáveis que estiverem na base serão exportadas no output, e não pode haver NA em nenhuma delas, por isso o comando "na.omit()" no código.

```{r base para correlação, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
Base_Correlacao <- Base_GECAT %>%
  dplyr::select("DA_w1", "DA.abs_w1", "SP", "POST", "Emerg", 
                "SIZE_w1", "ROA_w1", "BTM_w1","CFO_w1", "AU", "Law", "COV") %>%
  na.omit() 

view(Base_Correlacao)
```

## Forma simples entre duas variáveis

```{r}
cor(Base_Correlacao$SIZE_w1, Base_Correlacao$ROA_w1)
```

## Construir função

*O código a seguir é apenas para construir a função que calculará a correlação. NÃO pode ser alterado.*

```{r função de correlação, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
corstars <-function(x, method=c("pearson", "spearman"), removeTriangle=c("upper", "lower"),
                     result=c("none", "html", "latex")){
    #Compute correlation matrix
    require(Hmisc)
    x <- as.matrix(x)
    correlation_matrix<-rcorr(x, type=method[1])
    R <- correlation_matrix$r # Matrix of correlation coeficients
    p <- correlation_matrix$P # Matrix of p-value 
    
    ## Define notions for significance levels; spacing is important.
    mystars <- ifelse(p < .0001, "****", ifelse(p < .001, "*** ", ifelse(p < .01, "**  ", ifelse(p < .05, "*   ", "    "))))
    
    ## trunctuate the correlation matrix to two decimal
    R <- format(round(cbind(rep(-1.11, ncol(x)), R), 2))[,-1]
    
    ## build a new matrix that includes the correlations with their apropriate stars
    Rnew <- matrix(paste(R, mystars, sep=""), ncol=ncol(x))
    diag(Rnew) <- paste(diag(R), " ", sep="")
    rownames(Rnew) <- colnames(x)
    colnames(Rnew) <- paste(colnames(x), "", sep="")
    
    ## remove upper triangle of correlation matrix
    if(removeTriangle[1]=="upper"){
      Rnew <- as.matrix(Rnew)
      Rnew[upper.tri(Rnew, diag = TRUE)] <- ""
      Rnew <- as.data.frame(Rnew)
    }
    
    ## remove lower triangle of correlation matrix
    else if(removeTriangle[1]=="lower"){
      Rnew <- as.matrix(Rnew)
      Rnew[lower.tri(Rnew, diag = TRUE)] <- ""
      Rnew <- as.data.frame(Rnew)
    }
    
    ## remove last column and return the correlation matrix
    Rnew <- cbind(Rnew[1:length(Rnew)-1])
    if (result[1]=="none") return(Rnew)
    else{
      if(result[1]=="html") print(xtable(Rnew), type="html")
      else print(xtable(Rnew), type="latex") 
    }
} 
```

## Gerar tabelas

*Para análises com **variáveis categóricas** Spearman é o mais recomendado*

```{r gerar output excel, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
#pacotes usados: Hmisc, xtable

Matriz_Correlacao_Pearson <- corstars(Base_Correlacao,
                                      method = c("pearson"),
                                      removeTriangle = c("upper"),
                                      result = c("none", "html", "latex"))

Matriz_Correlacao_Spearman <- corstars(Base_Correlacao,
                                       method = c("spearman"),
                                       removeTriangle = c("lower"),
                                       result = c("none", "html", "latex"))


library(writexl)
write_xlsx(Matriz_Correlacao_Pearson,"Corr_pearson.xlsx")
write_xlsx(Matriz_Correlacao_Spearman, "Corr_spearman.xlsx")
```

# Regressão

## plm

### Estimando os modelos

*Variáveis Categóricas precisam da inclusão do "factor(VARIÁVEL)"*

**Pooled**

```{r pooled, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
mod_ACCY_Pooled <- plm(formula = ACCY_w1 ~ factor(POST)*factor(Emerg) + SIZE_w1 + ROA_w1 + Estim + factor(LCS),
                   data = Base_GECAT,
                   index = c ("Ticker","Year","Country"),
                   model = "pooling")

summary(mod_ACCY_Pooled)
```

**Efeitos Aleatórios**

```{r ea, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
mod_ACCY_EA <- plm(formula = ACCY_w1 ~ factor(POST)*factor(Emerg) + SIZE_w1 + ROA_w1 + Estim + factor(LCS),
                   data = Base_GECAT,
                   index = c ("Ticker","Year","Country"),
                   model = "random")

summary(mod_ACCY_EA)
```

**Efeitos Fixos**

O modelo "within" controla os efeitos fixos individuais;

O modelo "fd" estima interceptos específicos para cada grupo; e

O modelo "between" estima interceptos médios para cada grupo.

A escolha do modelo depende dos seus objetivos de pesquisa e da estrutura dos seus dados em painel.

within:

```{r within, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
mod_ACCY_EF_within <- plm(formula = ACCY_w1 ~ factor(POST)*factor(Emerg) + SIZE_w1 + ROA_w1 + Estim + factor(LCS),
                   data = Base_GECAT,
                   index = c ("Ticker","Year","Country"),
                   model = "within")

summary(mod_ACCY_EF_within)
```

fd:

```{r ft, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
mod_ACCY_EF_fd <- plm(formula = ACCY_w1 ~ factor(POST)*factor(Emerg) + SIZE_w1 + ROA_w1 + Estim + factor(LCS),
                   data = Base_GECAT,
                   index = c ("Ticker","Year","Country"),
                   model = "fd")

summary(mod_ACCY_EF_fd)
```

between:

```{r between, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
mod_ACCY_EF_between <- plm(formula = ACCY_w1 ~ factor(POST)*factor(Emerg) + SIZE_w1 + ROA_w1 + Estim + factor(LCS),
                   data = Base_GECAT,
                   index = c ("Ticker","Year","Country"),
                   model = "between")

summary(mod_ACCY_EF_between)
```

Comparando lado a lado os resultados

```{r ef lado a lado, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
export_summs(mod_ACCY_EF_within, mod_ACCY_EF_fd, mod_ACCY_EF_between)
```

Visualizando os resultados entre pooled, EA, e EF

```{r comparando modelos, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
export_summs(mod_ACCY_Pooled, mod_ACCY_EA, mod_ACCY_EF_within, mod_ACCY_EF_fd,  mod_ACCY_EF_between)
```

### Escolhendo o Melhor Modelo

#### Pooled x EF

Teste F de Chow \*A hipótese nula é de que há igualdade nos interceptos e nas inclinações para todos os indivíduos, caracterizando o modelo de dados agrupados (pooled). Se valor p\<0,05, o modelo de Efeitos Fixos é melhor do que o modelo Pooled.

```{r pooled x ef, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
pFtest(mod_ACCY_EF_between, mod_ACCY_Pooled) #EF
```

#### Pooled x EA

Teste de Breusch e Pagan \*A aceitação da hipótese nula implica que o modelo de dados agrupados (pooled) é preferível. Se valor p\<0,05, o modelo de Efeitos Aleatórios é superior ao modelo Pooled.

```{r pooled x ea, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
plmtest(mod_ACCY_Pooled, type="bp") #EA
```

#### EF x EA

Teste de Hausmann \*Se o teste rejeitar a hipótese nula, o modelo de Efeitos Fixos é o mais adequado. Se valor p\<0,05 o modelo de Efeitos Fixos foi considerado superior ao modelo de Efeitos Aleatórios.

```{r ef x ea, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
phtest(mod_ACCY_EF_between, mod_ACCY_EA) #EF
```

### Pressupostos

#### Distribuição dos Resíduos

Anderson-Darling normality test H0: Há distribuição normal dos resíduos - se p-value \< 0,05, rejeita-se a hipótese nula, logo não há distribuição normal dos resíduos.

```{r distribuição dos residuos, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
ad.test(mod_ACCY_EF_between$residuals)

hist(mod_ACCY_EF_between$residuals)
```

#### Multicolinearidade

Um VIF entre 5 e 10 indica alta correlação, o que pode ser problemático.

E se o VIF for acima de 10, você pode assumir que os coeficientes de regressão estão mal estimados devido à multicolinearidade.

```{r vif, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
car::vif(mod_ACCY_EF_between)
```

#### Heterocedasticidade

H0: homocedasticidade; p-value \> 0,05...

se p-valor \< 0,05 , rejeita-se a hipótese nula (homocedasticidade), logo há heterocedasticidade

H1: Heterocedasticidade; se p valor menor que 5%

Solução: rodar no lugar de mínimos quadrados generalizados, o mínimos quadrados ponderados -\> Solução: Método de mínimos quadrados ponderados que é um caso particular do método de mínimos quadrados generalizados (MQO), pode ser aplicado quando se diagnostica que a variância dos termos de erro depende da variável explicativa (Fávero & Belfiore).

```{r heterocedasticidade, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
bptest(mod_ACCY_EF_between)
```

#### Autocorrelação

H0: Ausência de autocorrelação entre os resíduos.

Com o p-value \< 0,05, podemos rejeitar a hipótese nula e concluir que existe autocorrelação entre os resíduos

```{r autocorrelação, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
pbgtest(mod_ACCY_EF_between)
```

## logit

### Montar base para logit só com as variáveis que quero

Todas as variáveis que estiverem na base serão exportadas no output, e não pode haver NA em nenhuma delas, por isso o comando "na.omit()" no código.

```{r base para logit, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
Base_Logit <- Base_GECAT %>%
  dplyr::select("Ticker", "Year", "Country", "DA_w1", "DA.abs_w1", "SP", "POST", "Emerg", 
                "SIZE_w1", "ROA_w1", "BTM_w1","CFO_w1", "AU", "Law", "COV", "LCS") %>%
  na.omit() 

view(Base_Correlacao)
```

### Estimando o modelo

```{r logit final, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
mod_SP_final <- pglm(factor(SP) ~ factor(POST)*factor(Emerg) + SIZE_w1 + ROA_w1 + 
                       BTM_w1 + CFO_w1 + factor(AU) + factor(COV),
                     family = binomial('logit'),
                     model = "pooling", data = Base_Logit,
                     index = c("Ticker", "Year", "Country"),
                     start=NULL,  method = "bfgs", print.level = 3, R = 5)

summary(mod_SP_final)
```

```{r logit, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
mod_SP <- glm(factor(SP) ~ factor(POST)*factor(Emerg) + SIZE_w1 + ROA_w1 +
                BTM_w1 + CFO_w1 + factor(AU) + factor(COV), 
              family = binomial('logit'), 
              data = Base_Logit)

summary(mod_SP)
```

```{r d, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
Base_Logit$probRestSP <- mod_SP$fitted.values
```

```{r conferir modelos logit, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
export_summs(mod_SP,mod_SP_final)
```

### Estimando a Razão de Chances

```{r razao de chance, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
logitor(factor(SP) ~ factor(POST)*factor(Emerg) + SIZE_w1 + ROA_w1 + 
          BTM_w1 + CFO_w1 + factor(AU) + factor(COV), 
        data = Base_Logit)
```

> **ROA_w1:**
>
> OR = 1.0190154
>
> Interpretação: Um aumento de 1 unidade em ROA_w1 está associado a um **aumento** de aproximadamente 1.9% nas chances de a empresa ter gerenciamento de resultados por pequenos ganhos (SP).
>
> **BTM_w1:**
>
> OR = 0.9627820
>
> Interpretação: Um aumento de 1 unidade em BTM_w1 está associado a uma **redução** de aproximadamente 3.7% nas chances de a empresa ter gerenciamento de resultados por pequenos ganhos (SP).
>
> **factor(AU)1:**
>
> OR = 1.9564506
>
> Interpretação: Estar na categoria 1 (BIG4) de AU está associado a um **aumento** de aproximadamente 95.6% nas chances de a empresa ter gerenciamento de resultados por pequenos ganhos (SP) *em relação à categoria de referência* (NBIG4)

#### Determinando o Intervalo de Confiança

\*A determinação do intervalo de confiança do modelo proposto é relevante para que seja analizada a estimativa do intervalo de predição do coeficiente da variável independente, a um nível de confiança de 95%. Desta forma, em 95% dos casos, o parâmetro dos coeficientes estará dentro deste intervalo.

```{r intervalo de confiança, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
exp(cbind(OR=coef(mod_SP_final), confint(mod_SP_final)))
```

### Qualidade do modelo

#### Curva ROC

```{r roc, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
require(pROC)

roc_SP=plot.roc(Base_Logit$SP,fitted(mod_SP))
```

```{r plot roc, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
plot(roc_SP,
     print.auc=TRUE, 
     auc.polygon=TRUE, 
     grud=c(0.1,0.2),
     grid.col=c("green","red"), 
     max.auc.polygon=TRUE, 
     auc.polygon.col="lightgreen", 
     print.thres=TRUE)
```

#### Matriz de confusão

Specificidade e Sensitividade

```{r matriz confusao, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
Base_Logit$pdata <- as.factor(ifelse(predict(mod_SP, newdata = Base_Logit, type = "response")>0.040,"1","0"))

confusionMatrix(Base_Logit$pdata, factor(Base_Logit$SP), positive="1")
```

# Exportar Resultados

Conhecendo os nomes das variáveis

```{r nomes para exportar regressão, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
export_summs(mod_ACCY_Pooled, mod_ACCY_EF_between, mod_ACCY_EA)
```

Exportando os modelos (pode ser apenas um ou mais outputs)

```{r exportar regressão, echo=TRUE, message=FALSE, warning=FALSE, , echo=TRUE, paged.print=FALSE}
export_summs(mod_ACCY_Pooled, mod_ACCY_EF_between, mod_ACCY_EA,
             scala = F,
             model.names = c("PAA_QA_Pooled","PAA_QA_EF","PAA_QA_EA"),
             coefs = c("Intercept" = "(Intercept)", 
                       "POST" = "factor(POST)1",
                       "Emerg" = "factor(Emerg)1",
                       "POST*Emerg" = "factor(POST)1:factor(Emerg)1", 
                       "SIZE" = "SIZE_w1", 
                       "ROA" = "ROA_w1", 
                       "Estim" = "Estim",
                       "Growth" = "factor(LCS)Growth", 
                       "Mature" = "factor(LCS)Mature", 
                       "Shake-out" = "factor(LCS)Turbulence", 
                       "Decline" = "factor(LCS)Decline"),
             error_pos = c("same", "below", "right"),
             bold_signif = 0.05,
             borders = 2,
             outer_borders = 2,
             statistics = c(N = "nobs.1", R2 = "r.squared",
                            adj.R2 = "adj.r.squared", p.value = "p.value", 
                            "GL" = "df", AIC = "AIC", "logLik" = "logLik",
                            "Pseudo R2" = "pseudo.r.squared"),
             scale = TRUE, robust = TRUE, digits = 3, vifs = TRUE,
             note         = "{stars}. Erros padrões robustos clusterizados",
             title = "Table X - Regression Models",
             to.file = "docx",
             file.name = "Output_Regressão.docx")
```

*\**Não importa a ordem que estavam as variáveis nos outputs, a ordem em que irão aparecer na tabela final é a mesma que for indicada no argumento "coefs".

\*Ao usar a função export_summs, "robust=true" indica que você deseja calcular estatísticas robustas para o modelo. Estatísticas robustas são usadas para fornecer inferências mais confiáveis quando as suposições de homocedasticidade (variância constante) e normalidade dos resíduos são violadas, o que pode ocorrer em modelos de regressão.
